MODULE fefmm_LIP_inner

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_init_inner,   &
             fefmm_free_inner,   &
             fefmm_inner,        &
             fefmm_trace

   REAL(REALK), ALLOCATABLE, SAVE :: Q(:)
   REAL(REALK), SAVE :: prefac
   LOGICAL, SAVE :: initialised = .FALSE.

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_init_inner(FEdim,LIPN,lambda)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN)  :: lambda
      INTEGER(INTK), INTENT(IN)  :: FEdim, LIPN

      ALLOCATE (Q(FEdim**3))

      SELECT CASE (LIPN)
         CASE (1)
            CALL fefmm_build_Q1(FEdim,lambda)
            prefac = half
         CASE (2)
            CALL fefmm_build_Q2(FEdim,lambda)
            prefac = one/3d0
         CASE (3)
            CALL fefmm_build_Q3(FEdim,lambda)
            prefac = 3d0/8d0
         CASE (4)
            CALL fefmm_build_Q4(FEdim,lambda)
            prefac = one/45d0
!         CASE (5)
!            CALL fefmm_build_Q5(FEdim,lambda)
!            prefac = 3d0/8d0
         CASE (6)
            CALL fefmm_build_Q6(FEdim,lambda)
            prefac = one/140d0
         CASE DEFAULT
            CALL fmm_quit('LIP NYI in fefmm_inner!')
      END SELECT

      prefac = (prefac*lambda)**3
      initialised = .TRUE.

   END SUBROUTINE fefmm_init_inner

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_free_inner

      IMPLICIT NONE
      DEALLOCATE (Q)
      initialised = .FALSE.

   END SUBROUTINE fefmm_free_inner

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_trace(func,trace)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN)  :: func(:)
      REAL(REALK),   INTENT(OUT) :: trace

      REAL(REALK) :: dummy(SIZE(Q))

      dummy(:) = one
      CALL fefmm_inner(func,dummy,trace)

   END SUBROUTINE fefmm_trace

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_inner(f1,f2,inner)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN)  :: f1(:), f2(:)
      REAL(REALK),   INTENT(OUT) :: inner

      INTEGER(INTK) :: k

      IF (.NOT. initialised) CALL fmm_quit('fefmm_inner not initialised!')

      inner = zero
      DO k = 1, SIZE(Q)
         inner = inner + Q(k)*f1(k)*f2(k)
      END DO
      inner = prefac * inner

   END SUBROUTINE fefmm_inner

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Q1(FEdim,lambda)
! Implemented for 1st order LIP functions

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim
      REAL(REALK),   INTENT(IN)  :: lambda

      INTEGER(INTK) :: i,j,k, ijk, hi,lo
      REAL(REALK)   :: faci,facj,fack

      !nd = 3*FEdim-2
      hi = 2*FEdim-1
      lo = FEdim

      ijk = 0
      DO i = lo, hi
         faci = two
         IF ( i == lo ) faci = one
         IF ( i == hi ) faci = one
         DO j = lo, hi
            facj= two
            IF ( j == lo ) facj = one
            IF ( j == hi ) facj = one
            DO k = lo, hi
               fack= two
               IF ( k == lo ) fack = one
               IF ( k == hi ) fack = one
               !ijk = k + (j-1)*nd + (i-1)*nd*nd
               ijk = ijk +1
               Q(ijk) = faci*facj*fack
            END DO
         END DO
      END DO

   END SUBROUTINE fefmm_build_Q1

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Q2(FEdim,lambda)
! Implemented for 2nd order LIP functions

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim
      REAL(REALK),   INTENT(IN)  :: lambda

      INTEGER(INTK) :: i,j,k, ijk, lo,hi,nd
      REAL(REALK)   :: faci,facj,fack

      !nd = 3*FEdim-2
      hi = 2*FEdim-1
      lo = FEdim

      ijk = 0
      DO i = lo, hi
         faci = get_fac(i)
         DO j = lo, hi
            facj = get_fac(j)
            DO k = lo, hi
               ijk = ijk +1
               !ijk = k + (j-1)*nd + (i-1)*nd*nd
               fack = get_fac(k)
               Q(ijk) = faci*facj*fack
            END DO
         END DO
      END DO

   CONTAINS

      FUNCTION get_fac(i)

         IMPLICIT NONE
         INTEGER(INTK), INTENT(IN) :: i
         REAL(REALK) :: get_fac

         get_fac = one
         IF ( i == lo ) RETURN   ! phi_1
         IF ( i == hi ) RETURN   ! phi_3

         IF ( MOD(i,2) == 0 ) THEN
            get_fac = 4.0d0    ! phi_2
         ELSE
            get_fac = two    ! phi_1 + phi_3
         END IF

      END FUNCTION get_fac

   END SUBROUTINE fefmm_build_Q2

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Q3(FEdim,lambda)
! Implemented for 3rd order LIP functions

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim
      REAL(REALK),   INTENT(IN)  :: lambda

      INTEGER(INTK) :: i,j,k, ijk, lo,hi,nd
      REAL(REALK)   :: faci,facj,fack

      nd = 3*FEdim-2
      hi = 2*FEdim-1
      lo = FEdim

      ijk = 0
      DO i = lo, hi
         faci = get_fac(i)
         DO j = lo, hi
            facj = get_fac(j)
            DO k = lo, hi
               ijk = ijk +1
               !ijk = k + (j-1)*nd + (i-1)*nd*nd
               fack = get_fac(k)
               Q(ijk) = faci*facj*fack
            END DO
         END DO
      END DO

   CONTAINS

      FUNCTION get_fac(i)

         IMPLICIT NONE
         INTEGER(INTK), INTENT(IN) :: i
         REAL(REALK) :: get_fac

         get_fac = one
         IF ( i == lo ) RETURN   ! phi_1
         IF ( i == hi ) RETURN   ! phi_4

         IF ( MOD((i-1),3) == 0 ) THEN
            get_fac = two    ! phi_1 + phi_4
         ELSE
            get_fac = 3.0d0    ! phi_2 or phi_3
         END IF

      END FUNCTION get_fac

   END SUBROUTINE fefmm_build_Q3

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Q4(FEdim,lambda)
! Implemented for 4th order LIP functions

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim
      REAL(REALK),   INTENT(IN)  :: lambda

      INTEGER(INTK) :: i,j,k, ijk, lo,hi,nd
      REAL(REALK)   :: faci,facj,fack

      nd = 3*FEdim-2
      hi = 2*FEdim-1
      lo = FEdim

      ijk = 0
      DO i = lo, hi
         faci = get_fac(i)
         DO j = lo, hi
            facj = get_fac(j)
            DO k = lo, hi
               ijk = ijk +1
               !ijk = k + (j-1)*nd + (i-1)*nd*nd
               fack = get_fac(k)
               Q(ijk) = faci*facj*fack
            END DO
         END DO
      END DO

   CONTAINS

      FUNCTION get_fac(i)

         IMPLICIT NONE
         INTEGER(INTK), INTENT(IN) :: i
         REAL(REALK) :: get_fac

         get_fac = 14d0
         IF ( i == lo ) RETURN   ! phi_1
         IF ( i == hi ) RETURN   ! phi_5

         IF ( MOD((i-1),4) == 0 ) THEN
            get_fac = 28d0      ! phi_1 + phi_5
         ELSE IF ( MOD((i),4) == 0 ) THEN
            get_fac = 64d0      ! phi_4
         ELSE IF ( MOD((i+1),4) == 0 ) THEN
            get_fac = 24d0      ! phi_3
         ELSE IF ( MOD((i+2),4) == 0 ) THEN
            get_fac = 64d0      ! phi_2
         END IF

      END FUNCTION get_fac

   END SUBROUTINE fefmm_build_Q4

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Q6(FEdim,lambda)
! Implemented for 6th order LIP functions

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim
      REAL(REALK),   INTENT(IN)  :: lambda

      INTEGER(INTK) :: i,j,k, ijk, lo,hi
      REAL(REALK)   :: faci,facj,fack

!      nd = 3*FEdim-2
      hi = 2*FEdim-1
      lo = FEdim

      ijk = 0
      DO i = lo, hi
         faci = get_fac(i)
         DO j = lo, hi
            facj = get_fac(j)
            DO k = lo, hi
               ijk = ijk +1
               !ijk = k + (j-1)*nd + (i-1)*nd*nd
               fack = get_fac(k)
               Q(ijk) = faci*facj*fack
            END DO
         END DO
      END DO

   CONTAINS

      FUNCTION get_fac(i)

         IMPLICIT NONE
         INTEGER(INTK), INTENT(IN) :: i
         REAL(REALK) :: get_fac

         get_fac = 41d0
         IF ( i == lo ) RETURN   ! phi_1
         IF ( i == hi ) RETURN   ! phi_7

         IF ( MOD((i-1),6) == 0 ) THEN
            get_fac = 82d0      ! phi_1 + phi_7
         ELSE IF ( MOD((i),6) == 0 ) THEN
            get_fac = 216d0     ! phi_6
         ELSE IF ( MOD((i+1),6) == 0 ) THEN
            get_fac = 27d0      ! phi_5
         ELSE IF ( MOD((i+2),6) == 0 ) THEN
            get_fac = 272d0     ! phi_4
         ELSE IF ( MOD((i+3),6) == 0 ) THEN
            get_fac = 27d0      ! phi_3
         ELSE IF ( MOD((i+4),6) == 0 ) THEN
            get_fac = 216d0     ! phi_2
         END IF

      END FUNCTION get_fac

   END SUBROUTINE fefmm_build_Q6

!-------------------------------------------------------------------------------

END MODULE fefmm_LIP_inner

